#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import os
import urllib
import numpy as np


def loadData():
    """ 载入数据 """
    # 下载数据文件
    dataFile = 'housing.data'
    if not os.path.isfile(dataFile):
        urllib.request.urlretrieve(
            'http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data',
            'housing.data')

    # 使用 numpy.loadtxt 载入数据
    data = np.loadtxt(dataFile)
    # 观察到数据中存在从某一位置开始全为0的列，做平滑处理
    data = data + 1e-9
    return data


def oneSamples(m=50):
    """ 预处理用于一元回归的样本数据 """
    data = loadData()
    # 取前m组数据
    x = data[0:m, 0]
    y = data[0:m, -1]

    # 排序
    index = np.argsort(x)
    x = x[index]
    y = y[index]

    # 标准化: (原始值 - 平均值) / 标准差
    x = (x - np.mean(x)) / np.std(x)
    y = (y - np.mean(y)) / np.std(y)

    # 数据转矩阵，不习惯二维数组
    X = np.matrix(x).transpose()
    Y = np.matrix(y).transpose()

    # 加入偏置单元 X <- [ 1, X ]
    X = np.hstack([np.ones((m, 1)), X])

    return (X, Y)


def leastSquare(X, Y):
    """ 最小二乘法计算参数矩阵 theta
    Parameters
    ----------
        X: 输入矩阵, m x (n + 1)
        Y: 输入矩阵, m x 1
    Returns
    -------
        theta: 参数矩阵, (n + 1) x 1
    """
    return np.linalg.pinv(X.transpose() * X) * X.transpose() * Y


def gradientDescent(theta,
                    X,
                    y,
                    alpha=1e-3,
                    lparameter=0,
                    maxIter=1000,
                    stopCost=1e-6):
    """ 
    梯度下降实现
    Parameters
    ----------
        X: 包含偏执单元的输入样本矩阵，大小为 m x (n + 1), m 为样本个数， n 为特征个数
        y: 输出样本矩阵, 大小为 m x 1
        alpha: 学习率，梯度下降的速率
        lparameter: 正规化参数
        maxIter: 最大迭代次数
        stopCost: 当损失小于此值时，提前终止迭代并返回结果
    Returns
    -------
        theta: 待求参数, (n + 1) x 1 的矩阵
        J: 损失
        counter: 迭代次数
    """
    counter = 0
    m = X.shape[0]
    J = np.zeros(maxIter)
    L = lparameter * np.eye(theta.shape[0])
    L[0][0] = 0

    while counter < maxIter:
        h = X * theta - y
        cost = 1 / (2 * m) * h.transpose() * h + \
            theta.transpose() * L * theta / 2
        J[counter] = cost
        counter += 1
        if cost < stopCost:
            break
        else:
            theta = theta - alpha / m * X.transpose() * (
                X * theta - y) + L * theta

    J = J[0:counter]

    return (theta, J, counter)
